---
title: "Experimental Evidence from Ecuador"
author: "Brett Bessen"
output: html_document
---

# Setup

```{r}

## Packages 
library(tidyverse)
library(conflicted)
library(psych)
library(here)
library(xtable)
library(marginaleffects)
library(texreg)
library(ggsci)
library(ggpubr)
library(wesanderson)
library(consort)
library(grid)


conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
conflict_prefer("mutate", "dplyr")

## Load data 

ecu_data <- read.csv(here("data/ecu_data.csv"))
```

## Additional Functions 

```{r}



#Function for adding labels
roundr <- function(x)
{ paste(format(round(x ,2), nsmall = 2)) }


## Plot themes --------------

brb_theme <- function() {
  theme_pubclean()  %+replace%
    theme( 
      axis.text.x = element_text(color = "black", size = 11),
      axis.text.y = element_text(color = "black", size = 11), 
      axis.title.y = element_text(color = "black", angle = 90),
      plot.title = element_text(hjust = 0.5), 
      strip.text.x = element_text(size = 11), 
      legend.position = "bottom")
} 


brb_theme_noleg <- function() {
  theme_pubclean()  %+replace%
    theme( 
      axis.text.x = element_text(color = "black", size = 11),
      axis.text.y = element_text(color = "black", size = 11), 
      axis.title.y = element_text(color = "black", angle = 90),
      plot.title = element_text(hjust = 0.5), 
      strip.text.x = element_text(size = 11), 
      legend.position = "none")
}


```



# Outcome Variables 

```{r}
## Recode DVs ----
ecu_data <- ecu_data %>% 
  mutate(should_close = (pop_post_1_1 - 1),               #should close council
         just_close = (pop_post_1_2 - 1),                 # justified to close
         close = (((pop_post_1_1 - 1) + (pop_post_1_2 - 1))/2), #close council
         dem_close = (exp_battery_4 - 1),                 #democratic to close council
         council_impeach = (exp_battery_1 - 1),           #council impeach mayor
         council_corrupt = (exp_battery_2 - 1),           # legislators are corrupt
         mayor_reps = (exp_battery_3 - 1))                #mayor represents people
```


# Main Independent Variables 

```{r}
## Recode round two vote  -----

ecu_data <- ecu_data %>%
  mutate(
    r2_vote = recode_factor(vote_r2, # round two vote
      `5` = "no vote",
      `1` = "rc",
      `2` = "creo",
      `4` = "null"),
    r2_v2 = recode_factor(vote_r2,
       `5` = "null/no vote",
       `4` = "null/no vote",
       `1` = "rc", 
       `2` = "creo")
  )

## Treatment: Mayors Arg ----
ecu_data$pop[ecu_data$pop==""] <- NA
ecu_data$pop <- factor(ecu_data$pop, 
                        levels = c("none", "corr", "full_pop"))


ecu_data$pop <- droplevels(ecu_data$pop)
table(ecu_data$pop)


## Mayor's Party Affiliation as a Factor ----
ecu_data$party[ecu_data$party==""] <- NA
ecu_data$party <- droplevels(factor(ecu_data$party))
table(ecu_data$party)

## Electoral support for Mayor's Party ---- 
ecu_data <- ecu_data %>% 
  mutate(align = case_when(r2_vote == "creo" & party == "creo"  ~ 1,
                           r2_vote == "rc"   & party == "rc"  ~ 1, 
                           r2_vote == "creo" & party == "rc" ~ 0, 
                           r2_vote == "rc"   & party == "creo" ~ 0,
                           r2_vote == "null" ~ 0, 
                           r2_vote == "no vote" ~ 0))


ecu_data$align <- factor(ecu_data$align)


# Non-voter/null-voter dummy:
ecu_data$non_null <- NA
ecu_data$non_null[ecu_data$r2_vote %in% c("rc", "creo")] <- 0
ecu_data$non_null[ecu_data$r2_vote %in% c("null", "no vote")] <- 1

ecu_data$non_or_null <- factor(ecu_data$non_null)

## Support Party, Opposition, or No Party ---
ecu_data <- ecu_data %>% 
  mutate(support = case_when(r2_vote == "creo" & party == "creo" ~ "Incumbent",
                             r2_vote == "rc"   & party == "rc" ~ "Incumbent", 
                             r2_vote == "creo" & party == "rc" ~ "Opposition",
                             r2_vote == "rc" & party == "creo" ~ "Opposition",
                             r2_vote == "null" ~ "None", 
                             r2_vote == "no vote" ~ "None"))


ecu_data$support <- factor(ecu_data$support, 
                            levels = c("Opposition", "None", "Incumbent"))


ecu_data <- ecu_data %>% 
  mutate(respondent_support = recode_factor(r2_vote,
                                            "null" = "none",
                                            "creo" = "creo",
                                            "no vote" = "none",
                                            "rc" = "rc"))

```

# Other Variables 

```{r}

## Executive Aggrandizement Attitudes
ecu_data <- ecu_data %>%
  mutate(ea_close_court = (exec_aggrand_1 - 10),
         ea_close_cong = (exec_aggrand_2 - 10),
         ea_limit_move = (exec_aggrand_3 - 10), 
         ea_decree = (exec_aggrand_4 - 10),
         ea_cni = (exec_aggrand_5 - 10))


ecu_data %>% 
  summarize(across(starts_with("ea"),
                   ~ mean(.x, na.rm = TRUE)))

ecu_data %>% 
  select(starts_with("ea")) %>% 
  psych::alpha()

ecu_data <- ecu_data %>% 
  mutate(
    ea =  ((ea_close_cong + ea_close_court + ea_limit_move + ea_decree + ea_cni)/5)
    )


```


```{r}
# Dataframe without missings for analysis

dat <- ecu_data %>%
  drop_na(pop, align, non_null, close)

```


# Descriptives

```{r}

# Table: Party Support by Discourse Condition
xtable(table(dat$respondent_support, dat$pop)) 


```


# Effect of Frame

```{r}
frame <- lm(close ~ pop, data = dat)

frame_preds <- predictions(frame, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop")))

ggplot(frame_preds) +
  geom_pointrange(aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = pop, y = estimate, label = roundr(estimate)), 
            nudge_x = .09, check_overlap = T) +
  scale_color_jama() +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist")) +
  labs(x = "Discourse Condition", y = "Mean Support", title = "Mean Support for Closing the Legislature") +
  brb_theme_noleg()

ggsave(file = here("figures/exp_fig1.pdf"),  height = 6, width = 8)

# Contrasts 
frame_contrasts <- comparisons(frame, variables = list(pop = "pairwise")) %>%
  summary()

frame_contrasts %>% 
  select(comparison = contrast, estimate, std.error, 
         statistic, p.value, conf.low, conf.high) %>% 
  xtable(digits = 2) %>%
  print(include.rownames = F)
```

# Conditional on Party Support

```{r}

cond_frame1 <- lm(close ~ pop*align + non_null, data = dat)

cond_frame1_preds <- predictions(cond_frame1, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop"), 
                                     align = c("0", "1")))

ggplot(cond_frame1_preds, aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, 
                              group = align, color = align, shape = align)) +
  geom_pointrange(position = position_dodge(width = .1)) + 
  geom_line(linetype = "dashed", show.legend = F, position = position_dodge(width = .1)) +
  scale_shape_discrete(name = "Party Support", labels = c("Opposition", "Incumbent")) +
  scale_color_manual(name = "Party Support", labels = c("Opposition",  "Incumbent"), 
                     values= wes_palette("Darjeeling1", n = 3)) +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist"))  +
  labs(y = "Mean Support", x = "Discourse Condition", title = "Mean Support for Closing the Legislature") +
  brb_theme()


# Save 
ggsave(
  here("figures/exp_fig2.pdf"), 
  height = 6, 
  width = 8)


#Effect of incumbent support by discourse 
cond_frame1_contrasts1 <- comparisons(cond_frame1, variables = "align", by = "pop")

#Out to latex
cond_frame1_contrasts1 %>%
  select(contrast, pop, estimate, std.error, 
         statistic, p.value, conf.low, conf.high) %>% 
  xtable(digits = 2) %>% 
  print(include.rownames = F)


# Effect of discourse by incumbent support
cond_frame1_contrasts2 <- comparisons(cond_frame1, variables = list(pop = "pairwise"), by = "align")

cond_frame1_contrasts2 %>%
  select(contrast, align, estimate, std.error, 
         statistic, p.value, conf.low, conf.high) %>% 
  arrange(contrast, align) %>% 
  xtable(digits = 2) %>% 
  print(include.rownames = F)
```

# Models to Latex 

```{r}

mods1 <- list(frame, cond_frame1)

screenreg(mods1)

screenreg(mods1, file = here("tables/exp_regs.txt"))

texreg(mods1, file = here("tables/exp_regs.tex"))

```

# Conditional Effects: Supporter/Opposition/Unaligned


```{r}
cond_frame3 <- lm(close ~ pop*support, data = dat)

cond_frame3_preds <- predictions(cond_frame3, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop"), 
                                     support = c("Opposition", "None", "Incumbent")))

ggplot(cond_frame3_preds) +
  geom_pointrange(aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, 
                      color = support, shape = support), 
                  position = position_dodge(width = .3)) +
  scale_shape_discrete(name = "Party Support") +
  scale_color_manual(name = "Party Support", 
                     values= wes_palette("Darjeeling1", n = 3)) +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist"))  +
  labs(y = "Mean Support", x = "Discourse Condition", title = "Mean Support for Closing the Legislature") +
  brb_theme()



# Save 
ggsave(
  here("figures/exp_fig3.pdf"), 
  height = 6, 
  width = 8)

```

# Manipulation Check: Mayor Represents People

## Main Effect 

```{r}
# Mayor represents the people 
mayor_reps <- lm(mayor_reps ~ pop, data = dat)


mayor_reps_preds <- predictions(mayor_reps, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop")))

ggplot(mayor_reps_preds) +
  geom_pointrange(aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = pop, y = estimate, label = roundr(estimate)), 
            nudge_x = .09, check_overlap = T) +
  scale_color_jama() +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist")) +
  labs(x = "Discourse Condition", y = "Mean Agreement", title = "Mayor Represents the Common People") +
  brb_theme_noleg()


# Save 
ggsave(
  here("figures/exp_fig4.pdf"), 
  height = 6, 
  width = 8)



# Contrasts
mayor_contrasts <- comparisons(mayor_reps, variables = list(pop = "pairwise")) %>%
  summary()

mayor_contrasts %>% 
  select(comparison = contrast, estimate, std.error, 
         statistic, p.value, conf.low, conf.high) %>% 
  xtable(digits = 2) %>%
  print(include.rownames = F)

```

## Conditional Effect 

```{r}

# Model
mayor_reps_cond <- lm(mayor_reps ~ pop*support, data = dat)
                        
# Predictions
mayor_reps_preds <- predictions(mayor_reps_cond, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop"), 
                                     support = c("Opposition", "None", "Incumbent")))
# Plot
ggplot(mayor_reps_preds) +
  geom_pointrange(aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, 
                      color = support, shape = support), 
                  position = position_dodge(width = .3)) +
  scale_shape_discrete(name = "Party Support") +
  scale_color_manual(name = "Party Support", 
                     values= wes_palette("Darjeeling1", n = 3)) +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist"))  +
  labs(y = "Mean Agreement", x = "Discourse Condition", title = "Mayor Represents the Common People") +
  brb_theme()


# Save
ggsave(
  here("figures/exp_fig5.pdf"), 
  height = 6, 
  width = 8)

# Contrasts 
mayor_reps_contrasts2 <- comparisons(mayor_reps_cond, 
                                     variables = list(pop = "pairwise"), 
                                     by = "support")
mayor_reps_contrasts2

mayor_reps_contrasts2 %>%
  select(contrast, support, estimate, std.error, 
         statistic, p.value, conf.low, conf.high) %>% 
  arrange(contrast, support) %>% 
  print(include.rownames = F)


```


# Heterogenous Treatment Effects

```{r}

# Heterogenous effects by prior exectutive aggrandizment attitudes
het_mod <- lm(close ~ pop*align*ea + non_null, data = dat)

# Predictions
het_preds <- predictions(het_mod, newdata = 
                            datagrid(pop = c("none", "corr", "full_pop"), 
                                     align = c("0", "1"),
                                    ea = c(0, 2, 4)))
# Relabel for plotting
het_preds <- het_preds %>% 
  mutate(ea = case_match(ea,
    0 ~ "Min. EA Index",
    2 ~ "Median EA Index",
    4 ~ "Max EA Index"))

# Relevel Factor
het_preds$ea <- factor(het_preds$ea, 
                       levels = c("Min. EA Index",
                                  "Median EA Index",
                                  "Max EA Index"))

# Plot Predictions
ggplot(het_preds, aes(x = pop, y = estimate, ymin = conf.low, ymax = conf.high, 
                              group = align, color = align, shape = align)) +
  geom_pointrange(position = position_dodge(width = .1)) + 
  geom_line(linetype = "dashed", show.legend = F, position = position_dodge(width = .1)) +
  scale_shape_discrete(name = "Party Support", labels = c("Opposition", "Incumbent")) +
  scale_color_manual(name = "Party Support", labels = c("Opposition",  "Incumbent"), 
                     values= wes_palette("Darjeeling1", n = 3)) +
  scale_x_discrete(labels = c("none" = "Baseline", "corr" = "Anti-Elitist", "full_pop" = "Populist"))  +
  facet_grid(~ea) +
  labs(y = "Mean Support", x = "Discourse Condition", title = "Support for Closing the Legislature") +
  brb_theme()

# Save 
ggsave(
  here("figures/exp_fig6.pdf"), 
  height = 6, 
  width = 8)

```


# Simulated Power Analysis

```{r}

set.seed(18481848)

# Set parameters
n_sims <- 1000    #  simulations
sample_sizes <- c(3013, 5000, 10000, 20000)  #sample sizes  

# Data frame for results 
sig <- data.frame(matrix(nrow = n_sims, ncol = length(sample_sizes)))
colnames(sig) <- c("n_3013", "n_5000", "n_10000", "n_20000")

for (j in 1:length(sample_sizes)) {
  
  sample_size <- sample_sizes[j]

  # Start the simulations
  for (i in 1:n_sims) {
  
    # Simulate data
    pop <-  factor(sample(c("none", "corr", "full_pop"), sample_size, replace = T))
    pop <- relevel(pop, ref = "none")
    align <-  factor(sample(c("0", "1"), sample_size, replace = T, prob = c(.65, .35)))
    non_voter <- factor(ifelse(align=="1", "0", 
                             sample(c("0", "1"), sample_size, replace = T, 
                                    prob = c(.7, .3))))
  
    # simulate y to assume the effect of populism is 50% larger among supporters
    # the rest of the parameters are based on estimates from cond_frame 1
    y <- 1.34 + .25 * as.numeric(pop == "corr") + .29 * as.numeric(pop == "full_pop") +
       0.20 * as.numeric(align== "1") + 0.16 * as.numeric(non_voter == "1") + 
       0.03 * as.numeric(pop == "corr") * as.numeric(align== "1") +
       0.10 * as.numeric(pop == "full_pop") * as.numeric(align== "1") + 
       rnorm(sample_size, sd  = 1.16)

    # Estimate Model
    model <- lm(y ~ pop*align + non_voter)
  
    # Tidy model 
    tidied <- tidy(model)
  
    # Save p-value from pop x align
    sig[i, j] <- tidied$p.value[7]
  }
}

power <- function(x) {
  x <- ifelse(x < .10, 1, 0)
  return(mean(x))
}

sapply(sig, power)

```

# Information on Coding Criteria

```{r}
all_obs <-  nrow(ecu_data) #3823 obs
nrow(dat) # 3013

# Didn't reach treatment  
drop_outs <- sum(is.na(ecu_data$pop))

# Missing on R2 vote
miss_vote <- ecu_data %>% 
  filter(!is.na(pop)) %>% 
  summarise(drop_r2 = sum(is.na(r2_vote)))


# Obs by treatment assignment
ecu_data %>% 
  filter(!is.na(pop) & !is.na(r2_vote)) %>% 
  group_by(pop) %>%  
  summarise(n = n())

  
# Missing on the DV by group
ecu_data %>% 
  filter(!is.na(pop) & !is.na(r2_vote)) %>% 
  group_by(pop) %>% 
  summarise(missing_dv = sum(is.na(close)))


# Distribution of support 
ecu_data %>% 
  drop_na(pop, r2_vote, close) %>% 
  group_by(pop, support) %>%  
  summarise(n = n()) %>% 
  mutate(total = sum(n))




```

# Consort Plot

```{r}
options(txt_gp = gpar(cex = .75))


txt1_side <- "Excluded (n=736):\n\u2022 Dropped out before randomization (n=727)\n\u2022 Missing 2nd round vote (n=9)"
txt2_side <- c("Excluded (n=26):\n\u2022 Missing on DV (n=26)", 
               "Excluded (n=30):\n\u2022 Missing on DV (n=30)", 
               "Excluded (n=18):\n\u2022 Missing on DV (n=18)")

txt3_final <- c("Final Analysis (n=1005): \n\u2022 Incumbent (n=358)  \n\u2022 Opposition (n=333)  \n\u2022 Non/Null (n=314)",
                "Final Analysis (n=1000): \n\u2022 Incumbent (n=345)  \n\u2022 Opposition (n=343)  \n\u2022 Non/Null (n=312)",
                "Final Analysis (n=1008): \n\u2022 Incumbent (n=335)  \n\u2022 Opposition (n=359)  \n\u2022 Non/Null (N=314)")

# Create Consort Plot
consort_plot <- add_box(txt = "Survey Sample (n=3823)")  %>% 
  add_side_box(txt = txt1_side)  %>%  
  add_box(txt = "Randomized (n=3087)")  %>% 
  add_split(txt = c("Baseline (n=1031)", "Anti-Elitist (n=1030)", "Populist (n=1026)")) %>% 
  add_side_box(txt = txt2_side)  %>% 
  add_box(txt = txt3_final) 

plot(consort_plot)

# Save
ggsave(
  file = here("figures/si_consort.pdf"), 
  plot = build_grid(consort_plot),
  height = 6, 
  width = 8)


```













